Disordered driven lattice gases with boundary reservoirs and Langmuir kinetics 
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The asymmetric simple exclusion process with additional Langmuir kinetics, i.e. attachment and 
detachment in the bulk, is a paradigmatic model for intracellular transport. Here we study this model 
in the presence of randomly distributed inhomogeneities ('defects'). Using Monte Carlo simulations, 
we find a multitude of coexisting high- and low-density domains. The results are generic for one- 
dimensional driven diffusive systems with short-range interactions and can be understood in terms 
of a local extremal principle for the current profile. This principle is used to determine current 
profiles and phase diagrams as well as statistical properties of ensembles of defect samples. 



I. INTRODUCTION 

Despite several recent investigations p], Q, d, 0, d, @, 0, 
i, i, ED, [H d El El EE El the influence of sitewise 
disorder in driven lattice gases is not yet fully understood 

One focus of studies on the influence of disorder and 
inhomogeneities was the Asymmetric Simple Exclusion 
Process (ASEP), especially its totally asymmetric vari- 
ant (TASEP). This process is not only believed to capture 
the essentials of driven diffusive syste ms, but its homoge- 
neous version is exactly solvable [H, El, Hfll • The exact 
solution allows to determine the steady state properties 
analytically without approximations. These results can 
then be used as reference system to study the influence 
of disorder, inhomogeneities etc. 

Here we will study the competition between disorder, 
realized through randomly distributed hopping rates as- 
sociated to the sites in the TASEP, and Langmuir kinet- 
ics, i.e. attachment and detachment processes in the bulk. 
This is not only of theoretical interest due to the chal- 
lenges posed by a non-trivial current profile, but also of 
direct relevance for the description of intracellular trans- 
port. The model which we will study here has originally 
been proposed to describe motor-based transport along 
microtubules. Although the microtubules itself are ho- 
mogeneous, the presence of microtubule- associated pro- 
teins (MAP) [21| can create inhomogeneities which in- 
fluence the motion of the motors [221 ]. 

In comparison to the ASEP, the current profile in the 
presence of Langmuir kinetics is no longer constant. This 
requires a slightly different approach since now a "local" 
point of view becomes necessary. Our main interest will 
be in the (local) transport capacity defined in Sec. [Ill This 
important observable is now also a local variable and is 
of direct relevance for biological applications. 

This paper is organized as follows: In Sec. [Ill we de- 
fine the models that are considered here and review some 
relevant results. Sec. [Til] reports results for current and 
density profiles obtained by computer simulations. In 
Sec. HVI we develop a theoretical framework that helps us 
to understand the simulation results and the phase dia- 



gram. This theoretical approach is applied in Sec. [V] to 
compute the probability that a randomly chosen defect 
configuration exhibits phase separation. Finally, Sec. EI 
gives a summary and conclusions. 



II. MODEL AND DEFINITIONS 

We consider driven lattice gases with open boundary 
conditions and Langmuir kinetics (LK). To be more spe- 
cific we focus on the TASEP which is believed to be a 
paradigmatic example for this class of dynamic processes. 
Here different extensions considering LK have been pro- 
posed, e.g. by including the diffusion of detached particles 
[23l [24|. We will focus on a less detailed model variant, 
the TASEP /LK [H, HJ, which is a TASEP with addi- 
tional particle creation and annihilation in the bulk. The 
TASEP is defined on a lattice of L sites which are num- 
bered from j = 1 to j = L beginning at the left. Each site 
can be occupied by at most one particle. The motion of 
the particles from left to right is defined by (local) tran- 
sition rates between adjacent sites. The corresponding 
hopping rates pj describing the transitions of particles to 
their right neighbours are inhomogeneous. We will focus 
on a binary distribution with two possible values pj = p 
or pj = q at each site j where q < p. Sites with transition 
rate pj = q will be referred to as defect sites, while a site 
with transition rate pj = p is called a non- defect site. In 
the following we will call a stretch of I consecutive defect 
sites a bottleneck of size I. 

The boundaries of the system are connected to reser- 
voirs so that particles can enter at the left end (j = 1) 
and leave at the right end (j = L). The (fixed) densities 
po and pL+i of the reservoirs control the effective entry 
and exit rates, a := po and (3 := 1 — Pl+i, respectively. 

Langmuir kinetics is realized by creation and annihali- 
tion of particles in the bulk. This can be interpreted as 
particle exchange with a bulk reservoir and corresponds 
to attachment and detachment processes in the biological 
context. The corresponding rates will be considered to be 
homogeneous, i.e. independent of the position, through- 
out this paper [4l|. 
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For large system size the investigation is usually sim- 
plified by performing a continuum limit. Since crucial 
properties, like the bottleneck lengths in a disordered sys- 
tem, might depend on the system size we have to specify 
this limit more carefully. We define a weak continuum 
limit where terms of 0(1/ L) are neglected while terms of 
Oil/ In L) are kept, and a strong continuum limit where 
we even neglect terms of 0(1/ In L). In the following we 
restrict ourselves to systems where the local creation and 
annihilation rates u a and cod are rescaled with the system 
size, while the global rates Vt a := oj a L and Qd '•= ^dL 
are kept constant. Hence Q, a and fid are system parame- 
ters while co a and uod are adjusted to the system size. In 
particular in the (weak and strong) continuum limit, the 
local rates vanish: w a , uJd — > for L — > oo. 

In homogeneous regions of these systems there is 
a unique current- density relation (CDR) J(p), usually 
called fundamental diagram in the context of traffic flow, 
that unambiguously gives the current for a given particle 
density p = (tj) on any site [27[, where Tj — 0, 1 is the 
occupation number of site j. The CDR of the TASEP 
has a single maximum. Later, when we will also consider 
more general driven lattice gases, we will always assume 
that their CDR also has a single maximum. The maxi- 
mum is at the point pM and takes the value Jm = J(pm)- 
In this case for a given current J, two possible values for 
the density, the high density value Ph(J) > Pm and the 
low density value Pl(J) < Pm exist. 

For these systems, the non-conservation of particles 
can be expressed by a source term in the equation of 
continuity of the stationary state [42^ : 



s(p) 



(1) 



where Jj is the current through the bond between sites 
j and j + 1. The attachment of particles is assumed to 
be inhibited by particles occupying sites, so we assume 
s(p) to be a globally decreasing function. In fact one 
can construct models with attractive interactions where 
s(p) is an increasing function. However, those systems 
might exhibit non-ergodic behaviour [28[ that we do not 
consider here. Since uo a , uod —> in the continuum limit, 
we also have s(p) — > in this limit. Hence locally the 
current is almost constant for large systems and the CDR 
is the same as in the corresponding system without LK 

nam. 

The time evolution per time interval At of the 
TASEP/LK can be written in terms of transition rules: 
For 1 < j < L: 



Hopping : 
Attachment : 
Detachment 

for j = 1: 



10 — > 01 with probability pjAt 

— > 1 with probability co a At (2) 

1 — > with probability UdAt 



Hopping : 10 
Entry : - 

Detachment : 1 — 



-> 01 with probability p\ At 
1 with probability a At (3) 
with probability UJdAt 



and for j = L: 

Attachment 
Exit : 



— > 1 with probability o; a At 

1 — > with probability (3 At 



(4) 



Other transitions are prohibited. Here "0" represents 
empty and "1" occupied sites. We can write the time 
evolution of the density pj = (tj ) as 



= Pj-liTj-liW ~ Tj {t))) - Pj { Tj (t)(l 
+ W a (l ~ Pj(t)) -U) d Pj(t) 



Tj+l(t))) 

(5) 



in the bulk and 
dpi 



dt 



dpL 
dt 



(t) = - Pl ( Tl (t)(l-r 2 (t))) 

+ a(l - pi(t)) -UdPi(t) 

(t) = PL-l(TL-l(t)(l-T L (t))) 

- (3p L (t)+uj a (l-p L (t)) 



(6) 



(7) 



at the left and right boundary, respectively. The param- 
eters a, (3 correspond to the generic boundary rates de- 
fined before. The source term is s(p) = co a (l — p) — LOdp- 
We call the hopping rates pj which are site-dependent 
properties intrinsic parameters which in the following 
will be considered as fixed, p = 1 and q = 0.6, if not 
stated otherwise. In contrast to this we consider the ex- 
plicit dependence of the system properties on the exter- 
nal parameters a, /?, ^ a and f^- Other driven lattice 
gases of the class characterized above can be written in 
the same way, while the local parameters might depend 
on the states in the vicinity of the sites and additional 
correlations might occur. Nonetheless one can assume 
that the TASEP/LK is quite universal as a paradigmatic 
model [H . 

In this work we are especially interested in randomly 
distributed defect sites. Here the defect density (j), which 
is the probability that a given site is a defect site, serves 
as an additional system parameter. Hence, transition 
rates are distributed as 



Pj 



q with probability (j) 
p with probability 1 



(8) 



Defect distributions of this kind are called disordered [43| . 
The properties of such systems are not fully determined 
by the defect density 0, but also depend on the spatial 
distribution of the defects. Since these properties can 
vary from sample to sample even for fixed system pa- 
rameters, an investigation of ensembles of systems (e.g. 
disorder average) rather than single samples is an issue 
of physical relevance. 

In the following sections we will make use of the 
particle-hole- symmetry exhibited by the TASEP/LK 
which is invariant under the symmetry operation 



0, a 



A 3 



i-j, n a ^n d . (9) 
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However, the particle-hole-symmetry is not essential for 
the generic behaviour, but it allows to reduce the param- 
eter space that needs to be investigated. 

The TASEP/LK with one defect site was already in- 
vestigated numerically and analytically in [30]. Now we 
want to generalize these results to arbitrary defect sam- 
ples. Therefor we introduce a local quantity, the trans- 
port capacity J*, which is the site-dependent maximum 
current that can be achieved by tuning the external pa- 
rameters a, /?, ft a and Qd in the continuum limit [441 ]. 
This quantity will be discussed in detail in section HVl 



III. OBSERVATIONS BY COMPUTER 
SIMULATIONS 
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In this section we summarize some properties of the 
system that can be observed with Monte Carlo simula- 
tions. Therefor we compare quantities of the inhomo- 
geneous TASEP/LK with the homogeneous TASEP/LK 
and the TASEP with defects. For simulations we used 
random-sequential update with fast hopping probability 
p = 1. If not specified else, we fix q = 0.6 as slow hop- 
ping rate. The unit of time is just one timestep so that 
probabilities and rates have the same numerical value. 



A. Few defects/ vanishing fraction of defects 

Before we consider finite defect densities <p > we dis- 
cuss systems with a fixed number of defects in the con- 
tinuum limit (0 = 0). Figs. [l]-[3] display the dependence 
of the densities and the current on the position in the 
system. 

Fig. [U shows the density and current profiles of a 
TASEP/LK-system with five defects, a homogeneous 
TASEP/LK-system and a TASEP with five defects in the 
low density phase. The density profiles of inhomogeneous 
and homogeneous TASEP/LK-systems differ only in the 
occurrence of narrow density peaks at the defects, while 
globally the density profile is the same. The current pro- 
files of the homogeneous and inhomogeneous system are 
identical. In contrast, the density profile of the TASEP 
with defects at the same sites shows density peaks as well, 
but the current profile (and the density profile far from 
the boundaries) is flat. This is due to particle conserva- 
tion while the lateral influx of particles allows a spatial 
variation of the current profile in the TASEP/LK where 
particles are not conserved in the bulk. 

Fig. [2] shows the corresponding situation for low exit 
rate and high entry rate. Due to particle- hole-symmetry, 
the results are analogous to the previous case. Adopting 
the terminology of the homogeneous system, the inhomo- 
geneous TASEP/LK-system can be considered to be in a 
high and low density phase, respectively. 

Fig. [3] displays density profiles for a « (3. As in the case 
above, homogeneous and inhomogeneous TASEP/LK- 
systems exhibit the same density profiles, apart from 
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FIG. 1: (Colour online) Comparison of current and den- 
sity profiles for a = 0.1 and (3 = 0.9 (low density phase) 
in the TASEP with defects, homogeneous TASEP/LK and 
TASEP/LK with defects and tt a = O d = 0.1. 



the peaks. In this case we see a shock in the density 
profile which is characteristic for non-particle-conserving 
dynamics in the bulk and which cannot be observed in 
the particle conserving TASEP (except at a = (3). 

Increasing the entry rate a for fixed and large (3 one 
observes a queuing transition in Fig.[H At a critical entry 
rate a* the peak at the leftmost defect broadens, forming 
a high density region. This corresponds to phase separa- 
tion and is also observed in the inhomogeneous TASEP 
at critical boundary rates. In the TASEP, however, the 
high density regime always extends to the left boundary. 
In contrast, the inhomogeneous TASEP/LK-system ex- 
hibits a stationary shock separating the low and high 
density region. Numerical finite-size scaling in Fig. [5] 
shows that the shock is getting sharper with increasing 
system size. Thus the high density region extends over a 
finite fraction of the system, corresponding to phase sepa- 
ration. In contrast, the peaks diminish for larger systems 
indicating that they are just local phenomena. We can 
associate this phase separation with a phase transition at 
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FIG. 2: (Colour online) Comparison of current and den- 
sity profiles for a — 0.9 and (3 = 0.1 (high density phase) 
in the TASEP with defects, homogeneous TASEP/LK and 
TASEP/LK with defects. 



the critical parameter value a' . 

Increasing a further moves the shock position to the 
left. The density profile right of the defect where phase 
separation occurred does not change anymore by varying 
the entry rate. The same is true for the output current 
at the right boundary J out = J(L). At some value of a 
a second high density region starts to form. Thus in a 
system with many defects multiple shocks can occur asso- 
ciated with alternating domains of high and low density. 

Above a critical value a*, where a high density do- 
main extends to the left boundary, the density profile 
and the current in the system is independent of the en- 
try rate. Since this independence also holds for large /?, 
we call this a Meissner phase in analogy to supercon- 
ductors, where the magnetic field in the interior bulk is 
independent of exterior fields. This terminology was also 
used for the boundary independent phase in the homo- 
geneous TASEP/LK [26[. However, one has to note that 
while in the homogeneous system there are long-range 
boundary layers in the density profile which do depend 



FIG. 3: (Colour online) Comparison of current and density 
profiles for a — 0.1 and f3 — 0.15 (high density phase) in the 
homogeneous TASEP/LK and TASEP/LK with defects. 

on boundary rates, the Meissner phase in the disordered 
system only exhibits short-range boundary layers. The 
current profile in fact does not depend on the boundary 
rates, both in the homogeneous and the inhomogeneous 
system. 

Due to particle-hole-symmetry all considerations made 
in this section can be transferred to the high density 
phase by replacing a with (3. 



B. Finite fraction of defects and disordered systems 

If the density of defects <fi is finite and the number of 
defects is of order of the system size, even a local in- 
crease of the density in the vicinity of the defects has 
considerable impact on the average density due to the 
large number of defects. The effect can be observed in 
Fig. [6] where we have simulated disordered systems with 
small but finite defect density <p for small a and large 
(3. In contrast to systems with few defects, the current 
profile of the disordered system differs from the that of 
the homogeneous system. This is due to the change of 
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FIG. 4: (Colour online) Density profiles for increasing values 
of a and fixed (3 = 0.9. At a critical value a* a high density 
region at the most right defect occurs (phase separation) . For 
higher a multiple high density regions appear. 
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FIG. 6: (Colour online) Comparison of current and density 
profiles for a — 0.1 and /3 = 0.9 (low density phase) in the 
disordered TASEP/LK with defect density = 0.1 and ho- 
mogeneous TASEP/LK. 



FIG. 5: (Colour online) Density profiles for identical macro- 
scopic parameters Q a = 0.1, = 0.1, a = 0.35, f3 = 0.9 
but different system sizes L. The left boundary of the high 
density region (shock) becomes steeper with increasing sys- 
tem size, indicating a macroscopic regime. 

the density by defects, which leads to an altered influx 
of particles in the bulk by attachment /detachment. So 
the gradient of the current profile in the disordered sys- 
tem is different from the one in the homogeneous system 
and also from the system with few defects because in the 
latter the effect on the average density is negligible. 

Like in the TASEP/LK with few defects we observe 
multiple high and low density domains for large boundary 
rates, which is displayed in Fig. [7| In fact it is harder to 
distinguish macroscopic high and low density regimes in 
the disordered case because of the rapid changes of den- 
sity on a microscopic scale. We have to simulate rather 
large systems in order to identify a macroscopic high (low) 
density domain by inspection. In Sec. Elwe introduce a 



numerical method that can detect high and low density 
domains automatically. 

IV. THEORETICAL TREATMENT 

In this section we develop a theoretical framework for 
the observations made by Monte Carlo simulations. We 
expect that concepts developed in this section are generic 
for a larger class of disordered driven lattice gases that 
have a single maximum in the current-density relation 
and weak induced effective interactions between defects. 
The restriction "weak interaction" is discussed in detail 
in [ll|. In addition, we assume that the bulk influx term 
S(p) is decreasing with increasing density. 

First we summarize the properties that distinguish 
the inhomogeneous (disordered) TASEP/LK from the 
TASEP and homogeneous TASEP/LK, respectively. 

1. In the TASEP/LK the particle number is not con- 
served in the bulk. Therefore generically the cur- 
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FIG. 7: (Colour online) Density profile for a — 0.9 and (3 — 
0.9 in the disordered TASEP/LK with defect density = 0.2. 
One observes phase separation with alternating high and low 
density domains. The black line displays the density, averaged 
over 50 adjacent sites. 



rent profile is not flat and stationary shocks can 
occur in the bulk. For particle conserving systems 
these are not possible [27|, |3l( . 

2. In the homogeneous TASEP, the current is re- 
stricted by the upper bound = p(l—p) = 0.25 
(for hopping rate p = 1) due to the bulk exclu- 
sion. Already a single defect site with lower hop- 
ping rate q < p reduces this maximum stationary 
current [32[. In [30[ it was shown that also in the 
TASEP/LK a single defect site d restricts the cur- 
rent by a value J* := J™ ax < J™£ at this site, that 
cannot be exceeded by tuning external parameters. 
The quantity is exactly the local transport ca- 
pacity defined in Sec. [Til However, due to the spa- 
tially varying current, this effect is only local and 
the maximum value of the current J™ ax on sites i 
far away from the defect can be larger than . For 
completeness, we define J* = on non-defect 

sites z, so the transport capacity is peaked on a 
single site. If the current imposed by the bound- 
ary rates is larger than the transport capacity of 
a defect, phase separation occurs, exhibiting sta- 
tionary shocks. In the inhomogeneous TASEP no 
stationary shocks can occur in the bulk, thus the 
high density regime always fills the whole system 
left of the current limiting defect. 



3. In systems with only few defects the relation be- 
tween the average density and the current at a 
given site is the same as in the homogeneous sys- 
tem. Thus current profiles are almost the same 
(as long as the maximum current is not exceeded). 
In disordered systems with a finite fraction of de- 
fects, however, the current-density relation is not 
the same as in the homogeneous system and de- 



pends on q and the distribution of defects, since 
the large number of density peaks have influence 
on the source term s(p) in (pQ) on a macroscopic 
scale. Therefore the current profiles differ from the 
homogeneous case. 

In order to capture these properties, we follow the con- 
cept of [3(j by focusing on the current profiles J(x). 



The influence of defects: 
conditions 



additional initial 



Locally the current profiles are determined by the con- 
tinuity equation (pQ). Introducing the continuous variable 
x := iJErfj which is the relative position in the system, 
one can write J;_i = Ji - \^ + 0(1/ L 2 ). In the sta- 
tionary state the continuity equation (pQ) becomes 



dJ 

dx 



S(p) + 0(l/L) 



(10) 



where the global source term S(p) = L s(p) was in- 
troduced. In the TASEP/LK, for example, we have 
S(p) = Q a (l — p) — fldP- I n the continuum limit we ne- 
glect terms of 0(1/ L) so that ([TO]) becomes an ordinary 
first order differential equation in the continuous variable 
x. The system, however, has at least two initial condi- 
tions (e.g. the boundary conditions in the homogeneous 
case), thus it is over deter mined. Each initial condition 
at a point xq is associated with one solution of the differ- 
ential equation ([T0|) J XQ (x) for the current and p XQ (x) for 
density, respectively. We call the mathematical solutions 
to single initial conditions J Xo (x) an d p Xo ( x ) local cur- 
rent/density profiles. Physically these solutions are not 
necessarily realized. 

For the TASEP/LK with a single defect it was shown 
by Pierobon et al. [Io[ that the finite transport capacity 
at the defect site, corresponding to a local upper bound 
of the current, can be regarded as an additional condition 
on the current profile. They argued that the local solu- 
tion of ([TO]) with the initial condition J(xd) = J*(xd) be- 
comes relevant if the local current profiles of the bound- 
ary conditions exceed J* at the defect site. Here we 
want to justify this approach and generalize it to a larger 
class of driven lattice gases with many defects, including 
randomly disordered systems, that meet the restrictions 
noted earlier in this section. 

In it was shown that the maximum current in par- 
ticle conserving driven lattice gases with randomly dis- 
tributed defects but low defect density depends approx- 
imately only on the size of the longest bottleneck (Sin- 
gle Bottleneck Approximation, SB A). This fact, together 
with the observations made in [30], motivates the gener- 
alization of the transport capacity to driven lattice gases 
(including TASEP/LK) with many defects but low defect 
density, introducing an approximation similiar to SBA. 
We call it the locally independent bottleneck approxima- 
tion (LIBA): The transport capacity at a site x, J*(x), 
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is approximately equal to the maximum current that can 
be achieved by tuning the boundary rates in the corre- 
sponding system containing only one bottleneck at this 
site |45j. Thus J*(x) can be obtained by refering to a 
single-bottleneck system where all other defects (except 
the bottleneck at site x) have been removed. 

In systems without LK the current is spatially con- 
stant and cannot exceed the minimum of J*(x) which 
corresponds to the transport capacity of the longest bot- 
tleneck, since in single bottleneck systems the maximum 
current is equal to the local transport capacity J* (x) and 
decreases with / [0, [l(J [33[ . In this case the LIB A reduces 
to the SBA. 

The LIBA neglects the influence of other defects on 
the transport capacity at site x. Nonetheless, we claim 
that the influence of other defects on the transport ca- 
pacity can be considered as a perturbation in the same 
way as it is the case for the SBA in particle conserving 
systems [10]. Since the local attachment and detach- 
ment rates vanish in the continuum limit, the transport 
capacity of a bottleneck should be the same as in the cor- 
responding particle conserving system. Therefore J*(x) 
is independent of Q a and f^. For the TASEP without 
LK analytical results are available [TO, [33[ that can be 
used to obtain approximations for the transport capac- 
ity. Since the maximal current in thes e sy stams depends 
only on the bottleneck length l(x) [a |_10[| this holds also 
for the transport capacity. The concept of a local trans- 
port capacity is applicable if interactions of defects near 
a bottleneck are not to large and distances of defects are 
not too small (i.e. low defect density Q). 

Hence, the transport capacity J*(x) yields an upper 
bound for the current profile, 

J(x)<J*(x) for all x, (11) 

while the function J*(x) of course is not continuous. 
Since on non-defect sites (which correspond to bottle- 
necks of size I = 0) the transport capacity is J* = Jy£££, 
it is sufficient to check condition ([TTJ) for defect sites. 
Their number is finite in finite systems but can be infi- 
nite in the continuum limit (e.g. for disordered systems 
with finite defect density). 

The problem of condition ([TT]) is that it is given as an 
inequality and does not provide initial conditions for ([TP]) 
on the defect sites. We now want to show that ([TT]) is 
identically fullfilled by a set of initial conditions 

J(x) = J*(x) at defect sites x , (12) 

if one assumes additionally that the physical local solu- 
tion at x is selected by shock dynamics. 

First of all, if we assume the conditions ([T2]) we see 
that, in contrast to the boundary conditions of the system 
which are usually given by a fixed density, the initial 
condition imposed by a defect provides the possibility 
of two realizations of the local density profile. Given 
the initial condition J(xo) = J*(xo) at a point xo, only 
the current is a fixed initial condition while, due to the 



non-unique inversion of the current-density relation (one 
maximum!), there are two possible values for the density, 
pH and pL (with pn > Pl)-, leading to two possible local 
solutions of ([TO]) , a high density solution Jh(x) and a low 
density solution Jl(x): 

* PH ► Jh(x ~ X , J*) 

J* / (13) 

PL ► Jl(x ~ Xo, J*) 

Taking into account shock dynamics, a constraint on the 
selection of a physical solution is given by the collective 
velocity 

v c (x) = J'(p(x)) (14) 

where J(p) is the current density relation and the prime 
denotes the derivative with respect to p [3J]. A solution 
can only propagate away from the initial point if the di- 
rection of v c is pointing away from it, i.e. left of it only 
solutions with v c < can exist, while right of it solutions 
must have v c > 0. In a system with a single maximum at 
density p m in the CDR, j£ > for p < p m and j£ < 
for p > p m , thus left of an initial point, only the high den- 
sity solution Jh can be realized, while right of it only Jl 
can physically exist. This principle is displayed in Fig. El 
top. Hence, each initial condition at a point xo can have 
its own solutions. We denote these local solutions by 

J(x-xo,r) = l J T H ^- Xo ' J *}f VX ^ X ° . (15) 
v u ' J \Jl(x-xo,J ) for x > xo v 7 

Actually the dependence on J* can easily be obtained by 
a shift operation if two functions Jl(x) and Jh(x) with 
initial conditions Jl(0) = j£ and Jij(0) = where j£ 
and are arbitrary chosen values in the high- and low 
density branch of the CDR. If the range in both branches 
of the CDR includes J = 0, one can simply choose j£ = 
j£ = [47]. Since the ODE (JTUD is of first order and 
does not explicitly depend on x, the high and low density 
solutions unambigiously depend on p and are monotonic. 
Thus different local solutions Jl,h can only differ by a 
shift in the variable x. An arbitrary solution Jl,h(x — 
xo, J*) can be obtained by shifting Jl,h(x) by an amount 
xl,h(J*) so that the value of the shifted function at x = 
is equal to J*. The functions xl,h(J*) are just the 
inverse functions of the unique functions Jl,h(x). Then 
the local solutions at a point with initial condition J* are 
given as 

J(x — xo, J*) = J(x — xo — x(J*)) • (16) 

The functions Jl,h(x) and xl,h(J) can for example be 
obtained by numerical solution of ([TU]) with initial con- 
ditions j£ H . 

B. Selection of the global current profile 

The physically realized global current profile in the 
steady state is also determined by shock dynamics [53, 
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FIG. 8: Top: Local solutions in the vicinity of a point with 
an initial condition J*. Due to the non-unique inversion of 
the current-density relation, there are two possible solutions. 
Since for a physical solution the direction of the collective ve- 
locity must point away from this position, only solutions with 
maximal current are realized. 

Bottom: Intersection point of local solutions of the density 
profile. The constraint that only upward shocks can exist im- 
plies that only solutions with minimal current are physically 
realized. 

l34j ]. Shocks manifest themselves as discontinuities in the 
density profiles. If they are stationary they connect dif- 
ferent local steady state solutions of ([TO]) to form a global 
solution. The crucial quantity for this selection is the 
shock velocity 



that determines the propagation of a discontinuity in a 
(not necessarily stationary) density profile. Here J + (p + ) 
is the current (density) right of the shock and J_ (p_) is 
the current (density) left of the shock. In homogeneous 
driven lattice gases with a single maximum in the CDR 
only upward shocks with p + > p_ can exist (see for ex- 



ample [H, [35[ ) . In this was generalized to systems 
with particle creation and annihilation in the bulk, as 
long as the local creation and annihilation rates vanish 
in the continuum limit, i.e. s(p) — > for L — > 00. In 
this case, the CDR is the same as in the corresponding 
particle conserving system. 

In inhomogeneous systems there can also be "down- 
ward" discontinuities at the defect sites due to the im- 
posed maximum current. However, these discontinuities 
usually are not called "shocks" since their dynamics dif- 
fer. In contrast to shocks they are sharp also in finite 
systems, thus there are no fluctuations. Due to the lo- 
cal character of v s and v c we can state that away from 
defects, where locally the system is homogeneous, only 
upward shocks can exist. 

Since the source term s(p) of (pQ) vanishes in the contin- 
uum limit, shocks can only be stationary at intersection 
points of a high and a low density solution Jh(x) and 
Jl(x). So only at these intersection points a switch of 
the physical realized local solution can occur. Note that 
local solutions of the same kind Jl or Jh cannot intersect 
since the differential equation ([TO]) is of first order. Since 
which determines the slope of the current profile, 
is assumed to be a monotonically decreasing function in 
p, we have S(pn) < S(pl), hence the gradient of the 
high density solution Jh(x) is smaller than the one of 
the low density solution Jl(x). Therefore left of an in- 
tersection point, we have Jl(x) < Jh(x), while right of 
it Jh(x) < Jl(x). Since Jl is the physical solution left 
of a shock and Jh right of it, always the minimal local 
solution is the physical one (see Fig. [HI bottom). We de- 
fine the minimal envelope of all the local current profiles 
as the capacity field of the system 

C(x) := min { J(x - x' , J*(x'))} (18) 

x' 

with defects at the points x' . This function does not 
depend on the boundary rates. The capacity field is a 
generalization of the capacity introduced in [30]. Note 
that in general the capacity field is not identical with the 
local transport capacity J*(x) [48]. The local transport 
capacity can be viewed as the source or "charge" of the 
capacity field. In this view, the function Jl,h(x — xo), 
which generates all local current profiles via ([TO]) , can be 
called the "Green's function" of the capacity field. 

Additional conditions on the current profile are given 
by the boundary rates so that p(0) = a and p(l) = 1-/3. 
Of course the maximum current of the homogeneous sys- 
tem J£^£ remains an upper bound also in the inhomo- 
geneous system. The capacity field together with the 
boundary conditions can be used to express the physi- 
cally realized current profile as 

J(x) =mm[J a (x),Jp(x),C(x)] (19) 

This principle is the generalization of the extremal cur- 
rent principle for the homogeneous TASEP [36[. It pro- 
vides a tool to obtain the global current profile if it is 
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possible to obtain the local solutions of ([TO]) and the lo- 
cal maximum current J*(x). 

Indeed the global current profile given by ([19)) iden- 
tically fulfills the condition ([T2|) that the current must 
always be lower than the transport capacity. 

In Fig. [9] we compare computer simulations of a system 
with a few defects with results obtained by the minimal 
principle ([T9|) in order to illustrate some features of the 
TASEP/LK with defects. We chose high boundary rates, 
so that the resulting current profile is exactly the capacity 
field C(x). For the values Q = Q a = Qd — 0.2 analyti- 
cal results for the local current profiles in the continuum 
limit are available. Following [25|, [3l| we used the refer- 
ence functions Jl(x) = Qx — Q 2 x 2 and Jh = —Flx + fl 2 x 2 
that obey the initial condition Jl,h(0) = to reproduce 
the local solutions of ([T0|) . The transport capacity was 
obtained in LIBA by results of a TASEP with a single 
bottleneck. The first three bottlenecks are well separated 
by a large distance. Here we see that LIBA works quite 
well and the current profile is reproduced by the minimal 
principle quite accurately. We also find that at the po- 
sition of bottleneck 2, the actual current is less than the 
transport capacity since the local solution of defect 1 is 
less than J*(x2) [49]. For bottleneck 4 there are devia- 
tions to LIBA since bottleneck 5 which is quite close to 
bottleneck 4 (distance = 6 sites) perturbs the transport 
capacity by further decreasing it. Nonetheless also in this 
region the minimal principle works if one takes the real 
transport capacity [50[ instead of LIBA. 
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C. Local current profiles in the disordered 
TASEP/LK 

We now want to quantify our results by finding the 
local solutions of the differential equation ([TO]) and the 
continuity equation (pp), respectively. For a numerical 
evaluation of these equations we need the CDR J(p) and 
its inverse Pl,h(J)- 

If there are only few defects in the system we have seen 
that the CDR is the same as in the homogeneous system, 
as long as the current is below the maximum current 
J*, since the increase of the average density is negligible. 
Thus in the TASEP/LK with defects we can use the same 
CDR as in the homogeneous system: J(p) = p(l — p). 
Therefore the local solutions are the same as the ones of 
the homogeneous systems. 

The situation is different for a finite fraction of defects 
in the system. Then the average density is strongly in- 
fluenced by the dense distribution of defect peaks which 
leads to an altered current-density relation even in the 
non-plateau region [6] . We will give an approximation to 
calculate the current-density relation for small, but finite, 
defect density <ft <C 1 if it is not too close to the maximum 
current. For that purpose we virtually divide the system 
into homogeneous subsystems with fast hopping rate p, 
while the slow hopping bonds connect these subsystems 
[5lj . In first instance we neglect correlations on the defect 



FIG. 9: (Colour online) Comparison of MC results and semi- 
analytical results for the capacity field (= current profile for 
high boundary rates; here a — (3 = 0.9) by LIBA. Bottlenecks 
are at sites x% (first defect site) with size h: 
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Further details are given in the text. 



bonds. The subsystems have an average size 1/0. In this 
point of view, the peaks at the defects are the boundary 
layers of the homogeneous subsystems. Without losing 
generality, we can assume the system to be in the low 
density phase and observe the local solution of the right 
boundary where peaks are concave. This can be trans- 
fered to high density solutions by particle-hole symmetry 
operation. Since uj a ^d ~ 1/L we can neglect them for 
large systems when looking at a single subsystem, thus 
we can treat them as homogeneous TASEPs. In a large 
homogeneous TASEP in the low den sity phase, the den- 
sity is given by po = 1/2 — y/T/i — J in the bulk far from 

the boundary. We can write the mass m := Y^t=i Pi °f 
the system as m = Lpo + m p with m p being the mass of 
the boundary layer. m p thus corresponds to the mass of 
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a peak in the inhomogeneous system. 

We approximate that the mass of the peaks does not 
depend on distance of adjacent defects. Then we can 
write the average density as 



p(x) = po(J(x)) + </>m p (J(x)), 



(20) 



since <\> is the fraction of defect sites. Surprisingly 
this rather uncontrolled approximation is supported by 
Fig. [10] where we plotted the mass in a system with two 
defects in dependence on the distance of latter ones. 
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FIG. 10: Mass m of two density-peaks m = ^2f =1 (pi — ot) 
in the low density phase of the TASEP with two defects in 
dependence on the distance of the defects. One observes that 
the dependence is rather weak. 

In this approximation, the mass of the peaks can be 
calculated analytically, since due to the independence of 
distance we can take it as the mass of the boundary layer 
in a large homogeneous TASEP, where exact results are 
available for given current J [18]. The density at a site 
L — n is given by 



(r L _ n ) = JS n (J) + J n+l R n (1/(1 - p)) (21) 



with 



R n {x) 



1 - VI -4x 



2x 



n+1 

£ 



(j-l)(2n-j)! 
n\(n + l- j)\ 



j=2 

Thus the peak mass is 

= i( T L-n) ~ a(l ~ a)] , 



(22) 
(23) 

(24) 



while the sum is truncated once the terms are small 
enough. 

Eqs. ([20|) - ([23|) can be used to calculate the current J 
for a given density p in the low density phase (and in the 



high density phase by particle-hole symmetry) and vice 
versa: 



J(p) = (p - (j>m p )(l - p + (j>m p ) . 



(25) 



This relation can be used to obtain a local solution of the 
differential equation ([TDD f° r a gi yen initial condition Ji 
by iteration. In Fig. [11] we compared profiles obtained by 
this procedure with results from computer simulations. 
One observes an excellent agreement which holds if the 
current is not close to the transport capacity. Together 
with the minimal current principle ([T9]) the global current 
profile can be obtained. 

The corresponding density profile can be obtained by 
inverting the CDR with respect to its two branches. Re- 
gions with a high density solution of the current profile 
correspond to a high density domain with the density 
Ph{J{ x )) obtained by the inverted current density re- 
lation. Analogous to that low density domains exist in 
regions of low density solutions. 
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FIG. 11: (Colour online) Comparison between simulation and 
analytical results for the current profile in a disordered system 
with cj) = 0.2, oo a = 0.2, cod = 0.1 for entry rate a = 0.1 and 
exit rate f3 — 0.9. Since the current is less than the transport 
capacity throughout the system the profile corresponds to the 
local current profile of the boundary condition p(0) = a. We 
observe excellent agreement between numerical and analytical 
results. This agreement holds for low current. Deviations 
occur only if the current comes close to the transport capacity. 



D. Phase diagram of disordered systems 

We now want to investigate the phase diagram of in- 
homogeneous driven lattice gases. 

If one of the local boundary solutions J a (x) or Jp(x) is 
the minimum of all local solutions in the whole system, 
we have a low density phase (L) in the former case and 
a high density phase (H) in latter one and there are no 
shocks in the system. These phases have the same macro- 
scopic properties like in the corresponding homogeneous 
system. 
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If there are intersecting points of local solutions they 
manifest themselves as shocks in the density profile, sep- 
arating high and low density regions (phase separation) 
corresponding to the realized high and low density solu- 
tions of the current profile. Phase separation can also be 
observed in homogeneous systems with Langmuir kinet- 
ics like the TASEP/LK and the NOSC-model 0, [33, [11 . 
Here the local solutions of the boundaries J a and Jp can 
intersect leading to a single stationary shock in the den- 
sity profiles, separating a low density domain left of it and 
a high density region right of it. This is called the shock 
phase (S) [25[ which is preserved as long as the minimum 
local profiles are the boundary current profiles. However, 
this kind of phase separation differs from the phase sepa- 
ration induced by defects. While in the S-phase the bulk 
behaviour is still determined by the boundary conditions, 
phase separation due to the finite transport capacity of 
defects is accompanied by a region where the current is 
"screened" by the defect (s) and is independent of the 
boundary condition, i.e. dJ Q^ > = for all x inside this 
region. If the phase separation is due to the screening by 
defects we rather refer to a defect-induced phase separated 
phase (DPS). If both boundary profiles J a (x) and Jp(x) 
are larger than C(x) in the whole system, the complete 
system is screened. The current profile is completely de- 
termined by the defect distribution and identical to the 
capacity field C(x). As argued in Sec. [Ill we call this fully 
screened phase Meissner phase (M). 

Another possible scenario is that the current near the 
boundaries is only limited by the maximum current of 
the bulk, i.e. C = and we have a maximum current 
phase with long ranging boundary layers like in the ho- 
mogeneous TASEP. However in disordered systems with 
randomly disordered defects, distances of defects are mi- 
croscopic and the probability that C = vanishes in 
the continuum limit. 

We can characterize the phases by two quantities: 

1. The total length A# of high density regions. This 
is the sum of individual high density regions and 
corresponds to the total jam length in traffic models 

& 

2. The screening length £ [52], which is the size of 
the area where the current profile does not depend 
on the boundary conditions. This is exactly the 
region where the boundary independent capacity 
field C(x) < Ja,p(%) and the local boundary profiles 
are not the physically realized ones. 

In tableUthe behaviour of these quantities in the different 
phases is displayed. Indeed this can be used to define the 
phases. For £ = defects do not influence the current 
profile and the system is in one of the "pure" phases, 
L,H or S, determined by the boundary conditions. If < 
£ < 1 there is phase separation and a part of the system 
does not depend on the boundary conditions, the system 
is in the DPS-phase. For £ = 1 the complete system 
is screened and the current profile is solely determined 
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TABLE I: Values and properties of the characteristic order 
parameters £ and A in the different phases. These properties 
can be used to define phases. 

by the defect distribution and the system is in the M- 
phase. The "pure" phases L,H,S can be characterized 
by £ = and the vanishing of high density regions (L, 
A = 0), coexistence of high and low density regions (S, 
< A < 1), and a global high density region (H, A = 1). 

The transition from L or H to DPS is marked by a 
discontinuity in £, but it is continuous in A. Indeed due 
to the discrete distribution of defects, £ itself is discon- 
tinuous throughout the DPS-phase while A is not. In the 
M-phase both £ and A are constant, while £ = 1 and A 
takes a finite value Am that is determined by the fraction 
of high density regions in the capacity field C(x) which 
depends on the individual defect distribution. 
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FIG. 12: (Colour online) Illustration of some current profiles, 
including critical profiles. We see that the critical rates are 
related by the critical current profiles: a* = pz,( J pi (1)), while 
Pl(p) is the inverted (low density) CDR and Jp/(1) is the lo- 
cal right boundary solution for (3 = [3 f . An analog relation 
is valid for f3* . The bold lines are the local current profiles 
consistent with the initial conditions imposed by the defects, 
whose minimal envelope is the capacity field. The thin lines 
are the critical boundary profiles and the dashed line corre- 
sponds to phase separated boundary current profiles. 

We see that at most phase boundaries both quantities £ 
and A are non-analytic. At the transition from S to DPS 
though A is analytic, thus it cannot be characterized by 
A. Hence for theoretical investigations it appears to be 
more convenient to use £ to discriminate defect- and non- 
defect phases. In simulations it is easier to detect phase 
separation (see next section) and use the non-analytic 
behaviour of A to obtain critical points. Due to the an- 
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alytic behaviour between S- and DPS-phase, however, 
this approach is only applicable at L-DPS and H-DPS- 
transitions. The S-DPS transition has to be obtained by 
theoretical considerations. 

In particle-conserving systems with defects the DPS- 
and S-phases vanish since no stationary shocks are possi- 
ble. Here both £ and A are discontinuous at the transition 
to the M-phase. However, in these systems the Meissner 
phase usually is also called "phase separated phase" [l(j 
since no distinction between several phases with phase 
separation has to be made. 

A sketch of the a — (3 phase diagram of a disordered 
driven lattice gas with LK is displayed in Fig.[l3l Attach- 
ment and detachment rates are fixed, while here uod > uo a . 
L-,H- and even S-phase might vanish for large fl a ,d if 
J a= o(L) > J* (£5) at some point Xd for any boundary 
rate a or (3 so that phase separation with screening al- 
ready occurs for vanishing boundary density. The dashed 
lines mark the phases of the homogeneous system. These 
pure phases are overlayed by the DPS- and M- phase 
which are characterized by the critical boundary rates 
a',/?' and a*, (3*. a' and /?' mark the minimal boundary 
rates at which the respective local boundary profile in- 
tersects the capacity field, i.e J a ^ > C(x) for at least one 
point x, while at the rates a*, /?*, J a ^ > C everywhere, 
so that local boundary profiles cannot propagate into the 
bulk. In Fig. [12] we sketched some critical current pro- 
files to illustrate the critical parameters. In parameter 
regions where J a and Jp do not intersect, a' and f3' do 
not depend on each other as well as a* and (3* , hence the 
phase diagram has a simple structure with phase bound- 
aries parallel to the parameter axes. Though, as we can 
see in Fig. [12j a' and (3* do depend on each other since 
J ((3*) = «/«'(!)• The same relation is valid for (3' and 
a*. Inside the region of intersecting boundary profiles 
(the shock phase of the homogeneous system), the struc- 
ture is nontrivial. The phase transition betweeen S- and 
DPS-phase depends explicitely on the variation of the 
intersection points of boundary profiles and minimal de- 
fect profiles. Explicitely it is given by the condition that 
a triple points x t with J a (xt) — Jpfat) = C(x) exist. One 
special case for which this condition can be solved exactly 
is the disordered TASEP/LK for Q a = in the strong 
continuum limit, where terms of 0(1/ InL) are neglected 
and the defect density <p scaling to zero as <\> ~ 1/lnL. 
In this case the capacity field C is constant and the tran- 
sition line is just a diagonal straight line. The phase 
diagram in the strong continuum limit is derived in the 
Appendix and displayed in Fig. [TH Though this limit is 
not quite physical it can be used as a reference point to 
argue that for finite defect densities the S-phase is convex 
(see also the Appendix). 

If we go away from the strong continuum limit, C(x) is 
not a constant. The structure of C is not smooth as was 
argued in Sec. IIVB[ so is the transition line. In Fig. [13] 
we displayed a rather generic sketch of a phase diagram 
that incorporates these arguments. Phase diagrams of 
other driven lattice gases with the properties noted in 
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FIG. 13: Phase diagram of the disordered TASEP/LK for 
> Hd. The critical rates depend on each other as 
a* = pl(J^(1)), which is argued in the text. The transition 
line between S- and DPS-phase is not smooth in the weak 
continuum limit due to the unsmooth structure of the capac- 
ity field (bold line). In the strong continuum limit the DPS 
phase is concave (bold dashed line). The topology of other 
disordered driven lattice gases is expected to be the same. 



the introduction will have the same topology. 



V. EXPECTATION VALUES FOR PHASE 
TRANSITIONS 

Like in particle conserving systems, the properties of 
disordered driven lattice gases with Langmuir kinetics de- 
pend strongly on microscopic details of the defect sam- 
ple. Since we are interested in macroscopic properties 
that do not depend on microscopic defect distributions, 
we concentrate on probabilistic quantities of ensembles of 
systems. One quantity of interest is the expected fraction 
of systems that exhibit phase separation in an ensemble 
of systems with identical system parameters and defect 
density. In this section we derive a procedure to calculate 
this quantity based on analytical results obtained by the 
principles from the last section. 

In order to compare these results with Monte Carlo 
simulations we introduce virtual particles similiar to sec- 
ond class particles [40[ that indicate if phase separation 
occurs in the simulated system. These particles do not 
change the dynamics of the system. The predicted prob- 
ability for phase separation is then compared with the 
relative frequency of phase separation in a set of simula- 
tions. 
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FIG. 14: Phase diagram of the disordered TASEP/LK with 
Q a = ='■ £1 in the strong continuum limit. The bold 
line at the S-DPS-boundary is valid for <j) scaling as 1/ In L 
and dashed line (sketched) is valid for finite defect density. 
The critical rates are given by a* = (3* = (1 — y/1 — q)/2, 
a' — a* — same for /3, with = 0.1, q = 15/16. The phase 
boundaries of the S-phase are of second order. For Q > 1/2, 
L- and H-phase vanish. 



A. Automated detection of phase separation 



B. Analytical approach for phase separation 
probability 

We use the results from the last sections in order to 
derive a analytical approach that allows the determina- 
tion of the probability that for a given defect density <\> 
phase separation occurs. Again we consider ensembles of 
systems instead of a fixed configuration of defects. 

The condition that no phase separation occurs is 

Ja(x) < J*(xb) and Jp(x) < J*(x&) for all x\>. (26) 

The fact that only low density solutions can intersect 
high density solutions also implies that an increases of 
a leads to a left shift of phase boundaries (in the phase 
separated phase) to the left while a increase of j3 moves 
the phase boundaries to the right. This can be seen in 
Fig. SI 

Following the LIBA we assume that the transport ca- 
pacity at a position x approximately depends only on 
the length of the bottleneck at this point, thus J*(x) ~ 
J* (/(#)). In a system with binary disorder there are on 
average L(l — 0) bottlenecks and the probability that one 
specific bottleneck has length I is P(l) = (1 — <p)(j) 1 [~u| . 

The relation between bottleneck length and transport 
capacity J* (I) as well as its inverse relation /(J*) can be 
obtained by analytical considerations or numerical com- 
putations in single bottleneck systems. The probability 
that the current is below the transport capacity at a given 
position x is then 



We introduce so-called virtual particles (V-particles) 
to identify and distinguish high and low density regions. 
These particles do not follow the exclusion constraint, 
instead they can occupy all sites even if these are oc- 
cupied by particles. The dynamics of the V-particles is 
the following: At the beginning, a V-particles is put on 
each defect site. After each lattice update the V-particles 
are updated sequentially beginning at the left. Each V- 
particle hops to the right if there is a particle on its site, 
while it hops to its left adjacent site if it is residing on an 
empty site. The V-particle cannot hop over slow bonds, 
thus if it is on a defect site, it cannot hop to the right, 
while if it is on a site right of a defect site, it cannot 
hop to the left. Hence, at any time, there is exactly one 
V-particle between each pair of contiguous defect sites. 
If the average density between two defects is larger than 
1/2, the V-particle tends to move to the right, while for 
p < 1/2 it tends to move to the left. Thus, we can 
identify a high density region by a V-particle that is, on 
average, closer to the right defect. By computing the av- 
erage distance of a V-particle to the defect right of it we 
can identify if there is a high density region in its vicinity. 

Using this procedure we can run a large number of 
simulations and automatically identify whether high and 
low density regions coexist. In this way the relative fre- 
quency of phase separated systems and an estimate for 
the probability of phase separation can be determined. 



p(j < j*) = p(i < i(j)) = p (o = 1 - LZ(J(x))J . 

l'=0 

(27) 
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FIG. 15: (Colour online) Fraction of samples that exhibit 
phase separation in dependence of the attachment rate uj a 
for fixed a = 0.1, /3 = 0.9, cod = 0.3. The system size is 
L = 1000 and each data point is obtained by simulating 200 
random defect samples with same system parameters. This is 
compared with analytical results obtained by ([28)) . 
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The probability V that no phase separation occurs is 
equal to the probability that the current is below the 
transport capacity everywhere in the system: 

(N bn ) L 

v= n p(j(i)< j*a(i))=na-^ (j(i))j ) • (28) 

i=l i=l 

Here N^n is the number of bottlenecks from left to right, 
so J(i) is the current at bottleneck i counted from the 
left. Since on average there are (Nb n ) = L(l — 0) bot- 
tlenecks, we can determine J(i) recursively by rescaling 
eq. (pQ) by the factor 1/(1 — 0) to obtain 

J(i + 1) = J(i)+o; a (l-0)(l-p(i))-^(l-0)p(i), 
p(i) = p (J(i))+(t)m p . (29) 

This way the probability for phase separation, which ex- 
plicitely depends on the system size can be computed 
iteratively by ([28]), while analytical results for J (I) in 
the TASEP with a single bottleneck are available [10|. 
In comparison to Monte Carlo simulations, this com- 
putation can be made with little effort. In Fig. [15] we 
simulated ensembles of random defect samples for differ- 
ent parameter values. The fraction of samples exhibiting 
phase separation is determined by the method from sub- 
section |VA] and compared with results obtained by ([28]) . 
One observes a region with a quite steep increase of the 
probability. The analytical results fit the simulation re- 
sults quite nicely, although there is a small shift to larger 
values of uo a . 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have investigated the interplay be- 
tween Langmuir kinetics (particle creation and annihi- 
lation in the bulk) and disorder, realized through ran- 
domly distributed hopping rates, in driven lattice gases 
connected to boundary reservoirs. Although both fea- 
tures provide a mechanism for phase separation (shock 
formation), the underlying mechanisms and dynamics is 
different and might lead to a form of competition. 

Based on Monte Carlo simulations of the disordered 
TASEP/LK, a TASEP with Langmuir kinetics and site- 
disordered hopping rates, the main properties of such 
systems have been identified. Like in the disordered 
TASEP we observe narrow peaks in the vicinity of defect 
sites. Their width however vanishes in the continuum 
limit. For larger values of the boundary rates we observe 
defect-induced phase separation, where multiple macro- 
scopic high and low density regions with a multitude of 
shocks occur. 

These findings can be understood in terms of an ex- 
tremal principle. In contrast to the principle originally 
proposed for homogeneous systems [34], [36j it is a local 
principle for the current profile. This is a direct conse- 
quence of the interplay between Langmuir kinetics, which 
induces a site-dependence of the stationary current, and 



the randomly distributed inhomogeneities. In our ap- 
proach we assumed that defects locally induce a reduced 
transport capacity imposing an upper bound to the cur- 
rent. For weakly interacting systems this quantity ap- 
proximately depends only on the local distribution of de- 
fects, especially on the size of the bottleneck. In this ap- 
proximation (LIB A) we can obtain the transport capacity 
by refering to single bottleneck systems. The transport 
capacity provides additional initial conditions to the dif- 
ferential equation ([TO]) that gives the slope of the local 
current profile in the continuum limit, each of them rep- 
resenting a individual local solution. Shock dynamics im- 
pose additional conditions on the physical current profile. 
Hence, out of the multitude of solutions only the profile 
that locally minimizes all solutions is physically realized. 
The full current profile can be obtained by superposing 
the solutions of all single bottlenecks which are described 
in terms of the same "Green's function" Jl,h(%) defined 
in Sec. HVBl 

While in systems with only few defects local current 
profiles are almost identical to those of the homogeneous 
systems, they significantly differ in large systems with 
a finite fraction of defect sites. In the case of the disor- 
dered TASEP/LK local density profiles can be accurately 
reproduced by identifying the density peaks with bound- 
ary layers of small virtual subsystems where exact results 
are available. 

The minimal principle can be used to predict some 
features of the phase diagram. As was already ob- 
served for single defects in [30[, defects can generate 
screened regions where the influence of boundary con- 
ditions vanishes. We can distinguish the original non- 
screened phases which are also present in homogeneous 
systems, a partially screened phase exhibiting phase sep- 
aration and a fully screened phase where the influence 
of the boundary conditions vanishes completely. For the 
strong continuum limit where terms of 0(1/ In L) do not 
contribute, the minimal principle even allows the deter- 
mination of the exact phase diagram, while in for the 
weak continuum limit at least most qualitative aspects 
of the phase diagram remain accessible. 

The LIBA and the minimal principle can also be ap- 
plied together with a statistical approach to obtain an 
approximation for the probability that a randomly pro- 
duced disorder sample exhibits phase separation. 

Although the results have been derived and tested on 
the TASEP/LK we believe that are generic for a large 
class of driven lattice gases, at least if they are ergodic 
with short-ranged interactions and a single maximum 
in the current-density relation. In more general pro- 
cesses the lateral current s(p) takes the role of attach- 
ment/detachment processes. 
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APPENDIX: PHASE TRANSITION LINES IN 
THE STRONG CONTINUUM LIMIT 

Usually it is quite difficult to determine the transi- 
tion line between S- and DPS-phase. One special case 
where it is possible to solve that problem exactly is in 
the strong continuum limit in the disordered TASEP/LK 
for Q a = Qd =: Q. In addition the number of defects 
is infinite, while the defect density is scaled to zero as 
(j) = 0(1/ In L). The average length of the longest bot- 
tleneck in a system of size L scales as In L/ In [III [l3[ , 
so in the strong continuum limit there has to be an in- 
finitely large bottleneck with a local transport capacity 
Jm = Moreover, we can say that this is the case for 
a ny sm all interval of length e if e is scaling slower than 
y/T/L, corresponding to LyJ\jL = ^/L sites. The global 
capacity field therefore simply is the constant function 
C(x) = g/4. Since the defect density vanishes, the CDR 
is the same as in the homogeneous system as was shown in 
the last sections numerically and analytically. The local 
boundary current and density profiles will therefore be 
the same as in the homogeneous system. Now the prob- 
lem we have to solve is equivalent to finding the transition 
from S- to LMH-phase in the homogeneous TASEP/LK 
if the homogeneous maximum current J* = 1/4 is ex- 
changed by q/4 [25|, |3l|. In these works, the transition 
line was determined to be (3* (a) = Pl{J*) — Ft — a. In- 
serting J* = g/4, we obtain for the transition line 



$*(&) = 1/2- 



1 



(30) 



which is just a shift of the pha se transition line to the 
right by the term y(l — qj/i. The properties of the 
phases of course are different to the ones in the homo- 
geneous system as we have argued before (especially the 
absence of long ranged boundary layers). The phase di- 
agram is displayed in Fig. [TH We have to point out that 



in this limit, the transition is of second order, since £ is 
continuous. 

Nonetheless the vanishing of (j) in the continuum limit is 
not quite physical, so we try to obtain at least qualitative 
results for the S-DPS transition line for finite (j). In sec. 
IIIIBI and llV CI we have seen that a small but finite defect 
density <j) > leads to a flattening of the local density 
profiles due to a broadening of the density peaks, so that 
their slopes dp Q X H , which are positive for Q a = a < 
1/2, (3 < 1/2, are decreasing for higher current J. 

Assume the system is on the transition line between S 
and DPS, i.e. a triple point x t with J a (xt) = Jp(%t) — 
q/4 exists. A shift of both p a and pp by an infinitise- 
mal amount dx also shifts the triple point though it per- 
sists. In parameter space, this corresponds to a move- 
ment along the transition line, while the boundary values 
are changed by 



da 



dp a 
dx 

dp 
da 



x=0 



dx and dft — — 

ox 



dx (31) 



x=l 



dp(3 




dx 


x=l 






dx 


x=0 



(32) 



using the relations a = p a (0) and (3 = 1 — pp(l). Since 
the boundary current J a ,/3 is monotonously increasing 
with a and (3 for a, (3 < 1/2, the flattening of the density 
profiles leads to: 



For (3 > a : 
For a > (3 : 



dp/3 
dx 

dpp 
dx 



< 



> 



dpg 

dx 

dpg 

dx 



x=0 



df3_ 
da 
dp 
da 



> -1(33) 
< -1(34) 



along the transition line. This corresponds to a concave 
distortion of the DPS-phase as displayed in Fig. [TH 
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